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Abstract 

We develop an incremental algorithm to compute the Newton polytope of the resultant, 
aka resultant polytope, or its projection along a given direction. The resultant is fundamental 
in algebraic elimination and in implicitization of parametric hypersurfaces. Our algorithm ex- 
actly computes vertex- and halfspace-representations of the desired polytope using an oracle 
producing resultant vertices in a given direction. It is output-sensitive as it uses one oracle 
call per vertex. We overcome the bottleneck of determinantal predicates by hashing, thus 
accelerating execution from 18 to 100 times. We implement our algorithm using the exper- 
imental CGAL package triangulation. A variant of the algorithm computes successively 
tighter inner and outer approximations: when these polytopes have, respectively, 90% and 
105% of the true volume, runtime is reduced up to 25 times. Our method computes instances 
of 5-, 6- or 7-dimensional polytopes with 35K, 23K or 500 vertices, resp., within 2hr. Com- 
pared to tropical geometry software, ours is faster up to dimension 5 or 6, and competitive in 
higher dimensions. 

Keywords General Dimension, Convex Hull, Regular Triangulation, Secondary Polytope, 
Resultant, CGAL Implementation, Experimental Complexity. 

1 Introduction 

Given pointsets Aq, . . . , An C Z", we define the pointset 

n 

A:=\J{Ax{e,})cZ^", (1) 

1=0 

where eo, . . . , e„ form an affine basis of M": cq is the zero vector, ei = (0, . . . , 0, 1, 0, . . . , 0), i — 
1, . . . ,n. Clearly, \A\ — \Ao\ + ■ ■ ■ + |A„|, where | • | denotes cardinality. By Cayley's trick (Prop.[2| 
the regular tight mixed subdivisions of the Minkowski sum + ' ■ ■ + An are in bijection with the 
regular triangulations of A, which are the vertices of the secondary polytope T,(A). 

The Newton polytope of a polynomial is the convex hull of its support, i.e. the exponent vectors 
of monomials with nonzero coefficient. It subsumes the notion of degree for sparse multivariate 
polynomials by providing more information, unless the polynomial is completely dense. Given n+1 
polynomials in n variables, with fixed supports Ai and symbolic coefficients, their sparse (or toric) 
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resultant TZ is a polynomial in these coefficients which vanishes exactly when the polynomials have 
a common root (Def. [T]). The resultant is the most fundamental tool in elimination theory, and 
is instrumental in system solving; it is also important in changing representation of parametric 
hypersurfaces. 

The Newton polytope of the resultant N{TZ), or resultant polytope, is the object of our study, 
especially when some of the input coefficients are not symbolic, in which case we seek a projection 
of the resultant polytope. The lattice points in N{TZ) yield a superset of the support of TZ; 
this reduces implicitization [T31 |3D] and computation of TZ to sparse interpolation (Sect. [2]). The 
number of coefficients of the n + 1 polynomials ranges from 0{n) to 0{n'^d'^), where d bounds 
their total degree. In system solving and implicitization, one computes TZ when all but 0{n) of 
the coefficients are specialized to constants, hence the need for resultant polytope projections. 

The resultant polytope is a Minkowski summand of Yj{A). For its construction, we exploit an 
equivalence relation defined on the secondary vertices, where the classes are in bijection with the 
vertices of the resultant polytope. This yields an oracle producing a resultant vertex in a given 
direction, thus avoiding to compute which has much more vertices than N(TZ). Although 

there exist efficient software for S(^) [28], it is useless in computing resultant polytopes. For 
instance, in implicitizing parametric surfaces with < 100 input terms, we compute the Newton 
polytope of the equations in < Isec, which includes all common instances in geometric modeling, 
whereas S(y^) is intractable. 

Our main contribution is twofold. First, we design an output-sensitive algorithm for computing 
the Newton polytope of TZ, or of specializations of TZ. The algorithm computes both vertex (V) 
and halfspace (H) representations, which is important for the targeted applications. Its incremen- 
tal nature implies that we also obtain a triangulation of the polytope, which may be useful for 
enumerating its lattice points in subsequent applications. The complexity is proportional to the 
number of output vertices and facets; the overall cost is dominated by computing as many regular 
triangulations of A (Thm. [lO| . We work in the space of the projected N{TZ) and revert to the 
high-dimensional space of only if needed. Our algorithm readily extends to computing S(^), 
projections of T,{A) and, more generally, any polytope that can be efficiently described by a vertex 
oracle. A variant of our algorithm computes successively tighter inner and outer approximations: 
typically, these polytopes have, respectively, 90% and 105% of the true volume, while runtime is 
reduced up to 25 times. This may lead to an approximation algorithm. 

Second, we describe an efficient implementation based on CGAL [7J and the experimental 
package triangulation. Our method computes instances of 5-, 6- or 7-dimensional polytopes 
with 35K, 23K or 500 vertices, respectively, in < 2hr. Our code is faster up to dimensions 5 or 
6, and competitive in higher dimensions, to a method computing N{TZ) via tropical geometry 
and based on the Gf an library [21) . Moreover, our code in the critical step of computing convex 
hulls, uses triangulation which compared to state-of-the-art software Irs, cdd, and polymake, 
is the fastest together with polymake. We factor out repeated computation by reducing the bulk 
of our work to a sequence of determinants: this is often the case in high-dimensional geometric 
computing. Here, we exploit the nature of our problem to capture the similarities of the predicates, 
and hash the computed minors which are needed later, to speedup subsequent determinants. This 
is of independent interest and can be used to improve other related problems such as convex hull 
and triangulation computations. 

Let us review previous work. Sparse (or toric) elimination theory was introduced in [18j . They 
show that N{TZ), for two univariate polynomials with kg + l, ki + 1 monomials, has {^°^''^) vertices 
and, when both ki > 2, it has fcofci -I- 3 facets. In [2^1 Sec. 6] is proven that N{TZ) is 1-dimensional 
iff \Ai\ = 2, for all i, the only planar N{TZ) is the triangle, whereas the only 3-dimensional ones 
are the tetrahedron, the square-based pyramid, and the polytope of two univariate trinomials; 
we compute an instance of the latter (Fig. [2]jb)). Following [29l Thm. 6. 2], the 4-dimensional 
polytopes include the 4-simplex, some N{TZ) obtained by pairs of univariate polynomials, and 
those of 3 trinomials, which we can investigate with our code: the maximal such polytope we have 
computed has f-vector (22, 66, 66, 22) (Fig. ^c)). 

In [5n] they describe all Minkowski summands of S(yl). In [57] is defined an equivalence 
class over S(A) vertices having the same mixed cells. The classes map in a many-to-1 fashion to 
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Figure 1: Let /q := Cqo - Coia:ia;2, /i := Cio - CiiXiX^, /2 := C20 - C2ia;f. Left: Newton polygons 
and those of fully dense polynomials of same total degree (dashed). Right: There are 2 regular 
tight mixed subdivisions of the Minkowski sum giving the extreme monomials oiTZ — — CqqC^ ]^C2i + 

CoiCloC20- 

resultant vertices; our algorithm exploits a stronger equivalence relationship. Tropical geometry is 
a polyhedral analogue of algebraic geometry which gives alternative ways of recovering resultant 
polytopes [21] and Newton polytopes of implicit equations [30]. Sect. [5] discusses it and shows our 
approach is faster up to dimensions 5,6 and competitive in higher dimensions. 

We focus on sequences of determinantal predicates. For determinants, the record bit complexity 
is 0(n^-^^'') [53]. Methods exist for the sign of general determinants, e.g. [5]. We compared linear 
algebra libraries LinBox [TT] and Eigen [T^, which seem most suitable in dimension > 100 and 
medium-high dimensions, respectively, whereas CGAL provides the most efficient determinant 
computation for the dimensions to which we focus. 

The roadmap of the paper follows: Sect. [2] describes the combinatorics of resultants, and the 
following section presents our algorithm. Sect. [4] overcomes the bottleneck of Orientation pred- 
icates. Sect. [5] discusses the implementation, experiments, and comparison with other software. 
We conclude with future work. 

2 Resultant polytopes and their projections 

We introduce tools from combinatorial geometry [551 [3T] to describe resultants [3 HI]. Let vol(-) € 
N denote normalized Euclidean volume, and (M™)^ the linear m-dimcnsional functionals. 

Let A C M."^ he a, pointset whose convex hull is of dimension d. For any triangulation T of A, 
define vector (j)T G M'-^' with coordinate 

</'T(a) = X! a&A, (2) 

aeT:aea 

summing over all simplices ct of T having a as a vertex; T,(A) is the convex hull of (jjr for all 
triangulations T. Let A"' denote pointset A lifted to via a generic lifting function w in 

(MI-^I)^. Regular triangulations of A are obtained by projecting the upper (or lower) hull of A^ 
back to W'^. 

Proposition 1. [18) The vertices of'S(A) correspond to the regular triangulations of A, while its 
face lattice corresponds to the poset of regular polyhedral subdivisions of A, ordered by refinement. 
A lifting vector produces regular triangulation T (resp. a regular polyhedral subdivision of A) iff it 
lies in the normal cone of vertex (pT (resp. of the corresponding face) o/S(^). The dimension of 
Y.{A) IS \A\-d-l. 

Let Aq, . . . , An be subsets of Z", Pqj • ■ • i -Pn C M" their convex hulls, and P — Pq + ■ ■ ■ + Pn 
their Minkowski sum. A Minkowski (maximal) cell of P is any full-dimensional convex polytope 
B — X]r=o^*' "^here each Bi is a convex polytope with vertices in Ai. Minkowski cells B,B' — 
^"^g B[ intersect properly when BiDB'^ is a face of both and their Minkowski sum descriptions are 
compatible, i.e. coincide on the common face. A mixed subdivision of P is any family of Minkowski 
cells which partition P and intersect properly. A cell is i-mixed or Vi-mixed, if it is the Minkowski 
sum of n one-dimensional segments from Pj, j ^ i, and some vertex Vi £ Pi. 



Mixed subdivisions contain faces of all dimensions between and n, the maximum dimen- 
sion corresponding to cells. Every face of a mixed subdivision of P has a unique description as 
Minkowski sum oi Bi <Z Pi- A mixed subdivision is regular if it is obtained as the projection of 
the upper (or lower) hull of the Minkowski sum of lifted polytopes P^^ := {{pi,Wi{pi)) \ pi G Pi}, 
for lifting Wi : — M. If the lifting function w := {wq . . . , Wn) is sufficiently generic, then the 
mixed subdivision is tight, and ^i^Q dim Bi = dim^"^Qi?i, for every cell. Given Aq,. .. ,An and 
the affine basis {eg, . . . , e„} of K", we define the Cayley pointset A C Z^" as in equation ([l]). 

Proposition 2. [Cayley trick] |18 l There exist bijections between: the regular tight mixed sub- 
divisions of P and the regular triangulations of A; the tight mixed subdivisions of P and the 
triangulations of A; the mixed subdivisions of P and the polyhedral subdivisions of A. 

The family A^,. . . ,An C Z" is essential if they jointly afhnely span Z" and every subset of 
cardinality j,l < j < n, spans a space of dimension > j. It is straightforward to check this 
property algorithmically and, if it does not hold, to find an essential subset. In the sequel, the 
input Aq, . . . , An C Z" is supposed to be essential. 

Definition 1. Let j4o,...,A„ C Z" be an essential family, and fo,..., fn G C[xi, . . . , Xn] 
polynomials with these supports and symbolic coefficients Cij,i — 0,...,n, j = 1,..., \Ai\, i.e. 
fi = X^aeA '^y*^"- Their sparse (or toric) resultant is a polynomial in 'L[cij : i — 0,...,n, 
j — I, . . . ,\Ai\], defined up to sign, which vanishes iff fo = fi = ■ ■ ■ = fn = has a common root 
(C*)", C* =C\{0}. 

For n = 1, the resultant is named after Sylvester. For linear systems, it equals the determinant 
of the (n + 1) X (n+ 1) coefficient matrix. The discriminant of a polynomial F{xi, . . . , Xn) is given 
by the resultant of F, dF/dxi, . . . , dF/dxn. 

The Newton polytope N(TZ) of the resultant is a lattice polytope called resultant polytope: 
the resultant has |.4| ~ X]r=o variables, hence N{TZ) lies in RI'^I, though of smaller dimension 
(Prop.[4|. The monomials corresponding to vertices of N{TZ) are the extreme resultant monomials. 
For a sufficiently generic lifting w G (M'-^I)^, the w-extreme resultant monomial, whose exponent 
vector maximizes inner product with w, is recovered from the regular tight mixed subdivision S 
of P induced by w. Let T be the regular triangulation corresponding, via the Cayley trick, to S, 
and pt G N'^' the exponent of the i/j-extreme monomial. Then, 

pria)^ ^ vo1((t) G N, a <G A, (3) 

a -mixed 

aeT:ae<y 

where simplex a is a- mixed iff the corresponding cell is a-mixed in S. Note that, pt{cl) G N, since 
it is a sum of volumes of mixed cells a, each of them equal to the mixed volume of a set of lattice 
polytopes, the Minkowksi summands of a: 

vol(cr) = MV{crQ, an), a = ctq -\ h (t„, (7^ C Ai, 

where MV denotes the mixed volume function. Now, N{TZ) is the convex hull of all pr vectors [l8l 

Eg. 

There exists a many-to-1 surjection from vertices of S(.4) to those of N(TZ). One defines 
an equivalence relationship on all regular tight mixed subdivisions, where equivalent subdivisions 
yield the same vertex in N{Tl). Thus, equivalent vertices of T,{A) correspond to the same resultant 
vertex. Consider w G (M'^')^ lying in the union of outer-normal cones of equivalent vertices of 
Ti{A). They correspond to a resultant vertex whose outer-normal cone contains w; this defines a 
w-extremal resultant monomial. If w is non-generic, it specifies a sum of extremal monomials in 
7^, i.e. a face of N{n) (Fig. [2];a),(b)). 

We compute some orthogonal projection of N{TZ), denoted 77, in M™: 

TT : MI-^I ^ M" : N{n) ^ 77, m < \A\. 
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(c) 



Figure 2: (a) E polytope of two triangles (dark, light grey) and one segment = 
{(0,0), (1,2), (4,1)}, Ai = {(0,1), (1,0)}, A2 = {(0,0), (0,1), (2,0)}; vertices correspond to mixed 
subdivisions of the Minkowski sum AQ + A1+A2 and edges to flips between them (b) N{TZ), whose 
vertices correspond to the dashed classes of S. Bold edges of S, called cubical flips, map to edges 
of N{TZ) (c) 4-dimensional N{TZ) of 3 generic trinomials with f-vector (22, 66, 66, 22); figure made 
with polymake. 

By reindexing, this is the subspace of the first m coordinates, so 7r(p) — (pi, . . . , p„i). It is 
possible that none of the coefficients Cij is specialized, hence m = |^|, tt is trivial, and U = N{TZ). 
Assuming the specialized coefficients take sufficiently generic values, U is the Newton polytope of 
the corresponding specialization of TZ. The following is used for preprocessing. 

Lemma 3. [2H Lem.3.20] If aij € Ai corresponds to a specialized coefficient of fi, and lies in the 
convex hull of the other points in Ai corresponding to specialized coefficients, then removing Oij 
from Ai does not change the Newton polytope of the specialized resultant. 

Proposition 4. [18j N(TZ) is a Minkowski summand ofT,{A), of dimension \A\ — 2n — 1. 

Let us describe the 2n + l hypcrplanes in whose intersection lies N{TZ). For this, let M be the 
(2n + 1) X matrix whose columns are the points in the Ai, where each a G Ai is followed by 
the i-th unit vector in N"+^. Then, the inner product of any coordinate vector of N{TZ) with row 
i of M is: constant, for i — 1, . . . , n, and known, and depends on i, for i = n + 1, . . . , 2n + 1 [181 
Prop. 7. 1.11]. This implies that one obtains an isomorphic polytope when projecting N{TZ) along 
2n + l points in A which affinely span M^" ; this is possible because of the assumption of essential 
family. Having computed the projection, we obtain N(TZ) by computing the missing coordinates 
as the solution of a linear system: we write the aforementioned inner products as M[X V]'^ = C , 
where C is a known matrix and [X V^]^ is a transposed (2n+l) x u matrix, expressing the partition 
of the coordinates to unknown and known values, where u is the number of N{TZ) vertices. If the 
first 2n + 1 columns of M correspond to specialized coefficients, M — [Mi M2], where submatrix 
Ml is of dimension 2n + 1 and inver tible, hence X = M^^{C - M2B). 

We focus on 3 applications. First, we interpolate the resultant in all coefficients, thus illustrat- 
ing an alternative method for computing resultants. 

Example 1. Let /g — a2X^ + aix + ag, fi — bix^ + bo, with supports Aq — 2,1,0, = 1,0. 
Their (Sylvester) resultant is a polynomial in 02, ai, Og, &i, &o- Our algorithm computes its Newton 
polytope with vertices (0,2,0,1,1), (0,0,2,2,0), (2,0,0,0,2); it contains 4 points, corresponding 
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to 4 potential monomials Uibibo, a^bi, a2aobibo, ai^o- There is a parameterization of resultants 
\24^ , which yields: 02 = (2ii +^2)^3^4; Q,i = (—2^1 — 2t2)t^t4^, oq = ^2^4; ^1 = — ^i^a^S; ^0 = ^i^Si 
where the ti 's are parameters. We substitute these expressions to the monomials, evaluate at 4 
sufficiently random ti 's, and obtain a matrix whose kernel vector (1, 1, —2, 1) yields TZ = a^bibo + 

agb^ — 2a2aofei6o + ai^o- ^ 

Second, consider system solving by the rational univariate representation of roots [3]. Given 
/i, ■ • • , /n G . . . , Xn], define an overconstrained system by adding /o = uo + uixi + • ■ • 

with symbolic u^'s. Let coefficients Cij,i > 1, take specific values, and suppose the roots of 
/i = ••• = /„ = are isolated, denoted = (r^i, . . . , ri„). Then the u-resultant is TZu = 
a Yir (""0 + uirn + • • • + Unrin)"^\ a G C* , where rrii is the multiplicity of r^. Computing TZu is 
the bottleneck; our method computes (a superset of) N{TZu)- 

The last application comes from geometric modeling, where yi = fi{x), i = 0, x = 

(xi, . . . ^Xn) G C M", defines a parametric hypersurface. Many applications require the equiva- 
lent implicit representation -F(yi, . . . ,yn) — 0. This amounts to eliminating x, so it is crucial to 
compute the resultant when coefhcients are specialized except the y^'s. Our approach computes 
a polytope that contains the Newton polytope of F, thus reducing implicitization to interpo- 
lation |14j . In particular, we compute the polytope of surface equations within fsec, assuming 
< 100 terms in parametric polynomials, which includes all common instances in geometric mod- 
eling, whereas the corresponding S(.4) are intractable. 



3 Algorithms and complexity 

This section analyzes our exact and approximate algorithms for computing orthogonal projections 
of polytopes whose vertices are defined by an oracle. This oracle computes the vertex of the 
polytope which is extremal in a given direction w. If there are more than one vertices extremal in 
w it computes exactly one of these. Moreover, we define such an oracle for the vertices of orthogonal 
projections U of N{TZ) which results to algorithms for computing 77 without computing N{TZ). 
Finally, we analyze the asymptotic complexity of these algorithms. 

Given pointset V, reg_subdivision(y, w) computes the regular subdivision of its convex hull 
by projecting the upper hull of V lifted by w, and conv(V^) computes the H-representation of its 
convex hull. The oracle Vtx(^, w, n) computes a point in 7T, extremal in the direction w, by 
refining the output of reg_subdivision(y^, w), w — {w,0) G (M'-^I)^, into a regular triangulation T 
of A, then returning Tr{px). It is clear that, triangulation T constructed by Vtx, is regular and 
corresponds to some secondary vertex (px which maximizes the inner product with Hi. 

Lemma 5. All points computed by Vtx are vertices of U . 

Proof. Let v = tt{pt) = Vtx(^, w;, tt). We first prove that v lies on 977. Point pr of N{TZ) is a 
Minkowski summand of vertex (f>T of S(^) extremal wrt w, hence px is extremal wrt w. Since w 
is perpendicular to projection tt, pt projects to a point in dU. The same argument implies that 
every vertex (jjip, where T' is a triangulation refining the subdivision produced by w, corresponds 
to a resultant vertex px' s.t. tt{pt') lies on the same face of 77 with px, hence also lies on 977. 

Now we prove that u is a vertex of 77. Let w be such that the face / of N{TZ) extremal wrt 
w contains a vertex px which projects to relint(7r(/)), where relint(-) denotes relative interior. 
Then, if Vtx(^, w,7r) returns tt{pt), this is a point on 977 but not a vertex of 77. We resolve 
such degeneracies by adding an infinitesimal generic perturbation vector to w, thus obtaining Wp. 
Since the perturbation is arbitrarily small, uip shall be normal to a vertex of /, extremal wrt 
w, but projecting to a vertex of 7r(/). This equivalent to computing a triangulation refining the 
regular subdivision of A induced by w, whose normal cone's projection is full (to) dimensional. 
The perturbation can be implemented in Vtx, without affecting any other parts of the algorithm. 
In practice, our implementation does avoid degenerate cases. □ 

The initialization algorithm computes an inner approximation of 77 in both V- and H-represen- 
tations (denoted Q,Q^ , resp.), and triangulated. First, it calls Vtx{A,w,TT) for w £ W C (M™)^; 
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set W is either random or contains, say, vectors in the 2m coordinate directions. Then, it updates 
Q by adding Vtx(^, w, tt) and Vtx(^, — w, tt), where w is normal to hyperplane H C containing 
Q, as long as either of these points lies outside H. We stop when these points do no longer increase 
dim(Q). 

Lemma 6. The initialization algorithm computes Q C 77 s.t. dim(Q) = dim(7T). 

Incremental Alg. [T] computes both V- and H-representations of U and a triangulation of 77, 
given an inner approximation Q, of 77 computed at initialization. A hyperplane 77 is legal 
if it is a supporting hyperplane to a facet of 77, otherwise it is illegal. At every step of Alg. [l] 
we compute v = Vtx(^, w, tt) for a supporting hyperplane 77 of a facet of Q with normal w. If 
w ^ 77, it is a new vertex thus yielding a tighter inner approximation of 77 by inserting it to Q, 
i.e. Q C CH(Q U u ) C 77, where CII(-) denotes convex hull. This happens when the preimage 
TT^^{f ) C N{Tl) of the facet f of Q defined by 77, is not a Minkowski summand of a face of 
having normal w. Otherwise, there are two cases: either w G 77 and v G Q, thus the algorithm 
simply decides hyperplane 77 is legal, or w G 77 and v ^ Q, in which case the algorithm again 
decides 77 is legal but also inserts v to Q. 

The algorithm computes from Q, then iterates over the new hyperplanes to either compute 
new vertices or decide they are legal, until no increment is possible, which happens when all 
hyperplanes are legal. Alg. [l] ensures that each normal w to a hyperplane supporting a facet of 
Q is used only once, by storing all used w's in a set W. When a new normal w is created, the 
algorithm checks if w ^ W, then calls Vtx(^, w, tt) and updates W -(^ W U w. If w G W then 
the same or a parallel hyperplane has been checked in a previous step of the algorithm. It is 
straightforward that w can be safely ignored; Lem. [7| formalizes the latter case. 

Lemma 7. Let 77' be a hyperplane supporting a facet constructed by Alg. [7| and H =/= H' an illegal 
hyperplane at a previous step. If H' , H are parallel then 77' is legal. 

The next lemma formulates the termination criterion of our algorithm. 

Lemma 8. Let vertex v be computed by Vtx{A, w, tt), where w is normal to a supporting hyperplane 
77 of Q, then v ^ H ^ H is not a supporting hyperplane of 11 . 

Proof. Let v — tt^pt), where T is a triangulation refining subdivision S in Vtx. It is clear 
that, since v S 977 extremal wrt w, if w ^ 77 then 77 cannot be a supporting hyperplane of 77. 
Conversely, let v £ H. By the proof of Lem. [5j every other vertex 7r(p^) on the face of N{TZ) is 
extremal wrt w, hence lies on 77, thus 77 is a supporting hyperplane of 77. □ 

We now bound the complexity of our algorithm. Beneath-and-beyond, given a /c-dimensional 
polytope with I vertices, computes its H-representation and a triangulation in 0{k^lt^), where t 
is the number of full-dimensional faces (cells) [32]. Let |77|, |77-^| be the number of vertices and 
facets of 77. 

Lemma 9. Alg. computes at most |77| + |77^| regular triangulations of A. 

Proof. The steps of Alg. [l] increment Q. At every such step, and for each supporting hyperplane 77 
of Q with normal w, we compute one vertex of 77, by Lem.|5] If 77 is illegal, this vertex is unique 
because 77 separates the set of (already computed) vertices of Q from the set of vertices of 77 \ Q 
which are extremal wrt w, hence, an appropriate translate of 77 also separates the corresponding 
sets of vertices of S(-4) (Fig. [s]). This vertex is never computed again because it now belongs to 
Q. The number of regular triangulations of A yielding vertices is thus bounded by |77|. 

For a legal hyperplane of Q, we compute one vertex of 77 that confirms its legality; the regular 
triangulation of A yielding this vertex is accounted for by the legal hyperplane. Since every normal 
to a hyperplane of Q is used only once in Alg. [T] (by the earlier discussion concerning the set W 
of all used normals), the statement follows. □ 
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Algorithm 1: Computeil {Aq, . . . , A„, tt) 

Input : essential Aq, . . . , An C Z" processed by lem. [3j 

projection tt : MI-^I E™, 

H-, V-repres. Q",Q; triang. Tq of Q C TT. 
Output: H-, V-repres. Q^,Q; triang. Tq oi Q ~ U. 

[JoiAi X a) II Cayley trick 

'^illegal ^ 0; 

foreach H e do HtUegai ^ TiiUegai U {ff} 
while Hiiiegai 7^ do 

select H g HiUegal and Huiegal ^ UtUegal \ {H}; 

w is tiie outer normal vector of H; 
V -i— Vtx(^, w, tt); 
ii V ^ HnQ then 

<3tonp^conv(QU{v}); 

foreach {d — l)-/ace f ^ Tq with f C 9iT do 
|_ Tq ^ Tq U {faces of conv{f, v)}; 

foreach iT' G {Q^ \ Q^^^} do 

Uuiegai <- Uiiiegai \ {H'jl I H' sepa.ra.tes Q , V 

1 

foreach i/' e {Q^^^ \ Q^} do 

Uaiegai ^ 'Huiegal U {H') II new h/plane 

return Q,Q",Tq; 



Let the size of a triangulation be the number of its cells. Let denote the size of the 
largest triangulation of A computed by Vtx, and sn that of U computed by Alg. [ij In Vtx, the 
computation of a regular triangulation reduces to a convex hull, computed in 0(n^|^|s^); for 
Pt we compute Volume for all cells of T in O(s^n^). The overall complexity of Vtx becomes 
0(n^|^|s3i). Alg. [ijcalls, in every step, Vtx to find a point on 977 and insert it to Q, or conclude 
that a hyperplane is legal. By Lem. [ojit executes Vtx as many as |7J| + |7T^| times, in 0((|7T| + 
|7J^|)n^|^|s^), and computes the H-representation of TT in 0{m^\n\sjj). Now we have, \A\ < 
{2n + l)sj^ and as the input |^|,to, n grows large we can assume that |TT| ^ \A\ and thus sn 
dominates s^. Moreover, 577(771+ 1) > |T7^|. Now, let O(-) imply that polylogarithmic factors 
are ignored. 

Theorem 10. The time complexity of Alg. [7] to compute 77 C M™ is 0{m^\II\s'jj + (|77| + 
\n"\)n^\A\s'^j^), which becomes 0(|7T|s|y) when |7T| > \A\. 

This implies our algorithm is output sensitive. Its experimental performance confirms this 
property, see Sect.jSj 

We have proven that oracle Vtx (within our algorithm) has two important properties: 

1. Its output is a vertex of the target polytope (Lem. [5]) 

2. When the direction w is normal to an illegal facet, then the vertex computed by the oracle 
is computed once (Lem. |9]) 

The algorithm can easily be generalized to incrementally compute any polytope P if the oracle 
associated with the problem satisfies property ([T]); if it satisfies also property ([2]), then the compu- 



8 



X 

Figure 3: Lem. [9] each illegal hyperplane of Q with normal w, separates the already computed 
vertices of U (here equal to N{JZ)) from new ones, extremal wrt w. X is a polytope s.t. X-\-N{Tl) — 
EM). 

tation can be done in 0(|P| + |-P^|) oracle calls, where |P|, |P^| denotes the number of vertices 
and number of facets of P, respectively. For example, if the described oracle returns Tr{(j)T) instead 
of Tr{pT), it can be used to compute orthogonal projections of secondary polytopes. 

The algorithm readily yields an approximate variant: for each supporting hyperplane H, we use 
its normal w to compute v =Vtx(^, w, tt). Instead of computing a convex hull, now simply take 
the hyperplane parallel to H through v. The set of these hyperplanes defines a polytope Qo ^ H, 
i.e. an outer approximation of U. Thus, we have an approximation algorithm by stopping Alg. [T] 
when yo\{Q)/vol{Qo) achieves a user-defined threshold. Then, vol((5)/vol(7I) is bounded by the 
same threshold. This algorithm is up to 25 times faster than the deterministic version. Of course, 
vol((5) is available by our incremental convex hull algorithm. However, vol(Qo) is the critical step; 
we plan to examine algorithms that update (exactly or approximately) this volume. 

When all hyperplanes of Q are checked, knowledge of legal hyperplanes accelerates subsequent 
computations of , although it does not affect its worst-case complexity. Specifically, it allows 
us to avoid checking legal facets against new vertices. 

4 Hashing of Determinants 

This section discusses methods to avoid duplication of computations by exploiting the nature of 
the determinants appearing in the inner loop of our algorithm. Our algorithm computes many 
regular triangulations, which are typically dominated by the computation of determinants. Similar 
techniques can be used to improve determinantal predicates in incremental convex hull computa- 
tions [13] . 

Consider the 2n x \A\ matrix with the points of A as columns. Define P as the extension of 
this matrix by adding lifting values w as the last row. We use the Laplace (or cofactor) expansion 
along the last row for computing the determinant of the square submatrix formed by any 2n + 1 
columns of P; wlog these are the first 2n+ 1 columns oi, . . . , a2n+i- Let (ai, . . . , a2n+i) \ i be the 
vector (ai, . . . ,a2n+i) without its i-th element; ^'(oi,...,a2„+i)\i is the (2n) x (2n) matrix obtained 
from the 2n first elements of the columns whose indices are in (ai, . . . , a2„+i) \ i. 

Orientation is the sign of the determinant of Pl^°™ +2) ' constructed by columns Oi , . . . , 02^+2 
when we add 1 S K^"+^ as last row. Computing a regular subdivision is a long sequence of such 
predicates with different a^'s. We expand along the next-to-last row which contains the lifting 
values and compute the determinants \P{ai,....a2n+2)\{i} \ fo'^ * ^ {!' . . . , 2n + 2}. Another predicate 
is Volume, used by Vtx. It equals the determinant of Pj^"^ constructed by columns 

ai, . . . , a2n+i when we replace its last row by 1 € M^"+-'^. 

Our contribution consists in maintaining a hash table with the computed minors, which will 
be reused at subsequent steps of the algorithm. We store all minors of sizes between 2 and 2n. 
For Orientation, they are independent of w and once computed they are stored in the hash table. 




N(n) '■■ EU) 
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The main advantage of our scheme is that, for a new w, the only change in P are m (nonzero) 
coordinates in the last row, hence computing the new determinants can be done by reusing hashed 
minors. This also saves time from matrix constructions. 

Laplace expansion has complexity 0(n!), which decreases as long as minors do not need to be 
recomputed. 

Lemma 11. Using hashing of determinants, the complexity of the Orientation and Volume pred- 
icates is 0{n) and 0(1), respectively, if all minors have already been computed. 

Many determinant algorithms modify the input matrix; this makes necessary to create a new 
matrix and introduces a constant overhead on each minor computation. Computing with Laplace 
expansion, while hashing the minors of smaller size, performs better than state-of-the-art algo- 
rithms, in practice. Experiments in Sect. [5] show that our algorithm with hashed determinants 
outperforms the version without hash. For m = 3 and m = 4, we experimentally observed that 
the speedup factor is between 18 and 100; Fig. |4(b)| illustrates the second case. 

We looked for a hashing function that takes as input a vector of integers and returns an integer 
that minimizes collisions. We considered many different hash functions, including some variations 
of the well-known FNV hash [T5]. We obtained the best results with the implementation of Boost 
Hash ED], which shows fewer collisions than the other tested functions. 

The drawback is the amount of storage, which is in 0(n!). The hash table can be cleared at any 
moment to limit memory consumption, at the cost of dropping all previously computed minors. 
Finding a heuristics to clear the hash table according to the number of times each minor was used 
would decrease the memory consumption, while keeping running times low. In our experiments, we 
obtained a good tradeoff between efficiency and memory consumption by clearing the entire hash 
table when it contains 10^ minors. Last column of Table [l] shows that the memory consumption 
of our algorithm is related to \A\ and dim{n). 

It is possible to exploit the structure of the above (2n) x (2n) minor matrices. Let M be such 
a matrix, with columns corresponding to points of Ag, . . . , A„. Let / : {1, . . . , 2n} — {0, . . . n} 
return i, given the column index. After column permutations, we split M into A nxn submatrices 
A, B, D, /, where / is the identity and I? is a permutation matrix. This follows from the fact that 
the bottom half of every column in M has at most one 1, and the last n rows of M contain at 
least one 1 each, unless detM = 0, which is easily checked. Now, det Af — ±det(_B — AD), with 
AD constructed in 0{n). Hence, the computation of (2n) x (2rt) minors is asymptotically equal to 
computing an n x n determinant. This only decreases the constant within the asymptotic bound. 
A simple implementation of this idea is not faster than Laplace expansion in the dimensions that 
we currently focus. However, this idea should be valuable in higher dimensions. 

5 Implementation and Experiments 

We implemented Alg. [TJin C++ to compute 77; our code can be obtained from http : //respol . ' 
sourcef orge .net. All timings are on an Intel Core 15-2400 3.1GHz, with 6MB L2 cache and 8GB 
RAM, running 64-bit Debian GNU/Linux. We rely on CGAL, using mainly a preliminary version 
of package trismgulation under development, for both regular triangulations, as well as for the 
V- and H-representation of 77. 

We first compare 4 state-of-the-art exact convex hull packages, triangulation 4 imple- 
menting [S], and beneath-and-beyond (bb) in polymake |17) : double description implemented 
in cdd |16j . and Irs implementing reverse search [1], to compute 77 actually extending the work 
in [2] for the new class of polytopes 77. The first was shown to be faster in computing Delau- 
nay triangulations in < 6 dimensions [1]. The last 3 are run through polymake, where we have 
ignored the time to load the data. We test all packages in an offline version. We first compute 
the V-representation of 77 using our implementation and then we give this as an input to the 
convex hull packages that compute the H-representation of 77. Moreover, we test triangulation 
by inserting points in the order that Alg. [l] computes them, while improving the point location of 
these points since we know by the execution of Alg. fllone facet to be removed (online version). 
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TTi 


1 A\ 


#of7T 


time 


respol 


vertices 


respol 


tr/on 


tr/off 


bb 


cdd 


Irs 


memory 


3 


2490 


318 


85.03 


0.07 


0.10 


0.07 


1.20 


0.10 


37 


4 


27 


830 


15.92 


0.71 


1.08 


0.50 


26.85 


3.12 


46 


4 


37 


2852 


97.82 


2.85 


3.91 


2.29 


335.23 


39.41 


64 


5 


15 


510 


11.25 


2.31 


5.57 


1.22 


47.87 


6.65 


44 


5 


18 


2584 


102.46 


13.31 


34.25 


9.58 


2332.63 


215.22 


88 


5 


24 


35768 


4610.31 


238.76 


577.47 


339.05 


> Ihr 


> Ihr 


360 


6 


15 


985 


102.62 


20.51 


61.56 


28.22 


610.39 


146.83 


2868 


6 


19 


23066 


6556.42 


1191.80 


2754.30 


> Ihr 


> Ihr 


> Ihr 


6693 


7 


12 


249 


18.12 


7.55 


23.95 


4.99 


6.09 


11.95 


114 


7 


17 


500 


302.61 


267.01 


614.34 


603.12 


10495.14 


358.79 


5258 



Table 1: Total time (in seconds) and memory consumption (in megabytes) of our code (respol) 
and time comparison of online version of triangulation (tr/ on) and offline versions of all convex 
hull packages for computing the H-representation of 77. 



The experiments show that triangulation and bb are faster than Irs, which outperforms cdd. 
Furthermore, the online version of triaingulation is 2.5 times faster than its offline counterpart 
due to faster point location (Table [T]). 



input 


m 
1^1 


3 3 4 4 5 5 
200 490 20 30 17 20 


approxim. 
algorithm 


# of Q vertices 
vol(Q)/vol(i7) 
vol(Qo)/vol(i7) 
time 


15 11 63 121 > lOhr > lOhr 
0.96 0.95 0.93 0.94 > lOhr > lOhr 
1.02 1.03 1.04 1.03 > lOhr > lOhr 
0.15 0.22 0.37 1.42 > lOhr > lOhr 


uniformly 
random 


\Q\ 

random vectors 
vol(Q)/vol(7T) 
vol(go)/vol(7T) 
time 


34 45 123 207 228 257 
606 576 613 646 977 924 
0.93 0.99 0.94 0.90 0.90 0.90 
1.05 1.01 1.02 1.02 1.03 > lOhr 
5.61 12.78 1.1 4.73 8.41 16.9 


random 
variant 


# of Q vertices 
vol(Q)/vol(i7) 
vol(Qo)/vol(i7) 
time 


26 21 102 380 341 544 
0.92 0.90 0.80 0.92 0.80 0.81 
1.02 1.04 1.14 1.03 1.08 1.10 
0.16 0.31 1.54 23.66 59.87 211.50 


exact 
alg. 


# of 77 vertices 
time 


98 133 416 1296 1674 5093 
2.03 5.87 3.72 25.97 51.54 239.96 



Table 2: Typical results on experiments computing Q, using the approximation algorithm and 
the random vectors procedures (the uniform and its variant); we stop the approximation algorithm 

when yQY{Q^) ^ results on random procedures is the average values over 10 independent 
experiments; "> lOhr" means computation of vo\{Qo) takes > lOhr on polymake. 

An advantage of triangulation is that it maintains a polytope whose boundary and interior 
are triangulated. This is useful when the oracle Vtx(^, w, tt) needs to refine the regular subdivi- 
sion of A which is obtained by projecting the upper hull of the lifted pointset (Sect. [3]). In 
fact this refinement is attained by a placing triangulation, i.e., by computing the projection of the 
upper hull of the placing triangulation of A"' . This is implemented in two steps: 

Step 1. compute the placing triangulation of the last \A\ — m points in the order they appear in 
A^ (they all have height zero), 

Step 2. place the first m points of A^ in To in the order they appear in A. 
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Figure 4: (a) Implicitization and M-resultaiits for 71 — 2,171 — 3 (b) Comparison of respol 
(hashing and not hashing determinants) and Gf an (traversing tropical resultants and computing 
normal fan from stable intersection) for m = 4 (c) Performance of Alg [l] for m = 3, 4, 5 as a 
function of input (d) Performance of Alg [l] as a function of its output; y-axis in (b),(c),(d) are in 
logarithmic scale. 



Step 1 is performed only once at the beginning of the algorithm, whereas Step 2 is performed every 
time we check a new w. The order of placing the points in Step 2 only matters if w is not generic; 
otherwise, w already produces a triangulation of the m points, so any placing order results in this 
triangulation. 

We can formulate this 2-stcp construction using a single lifting. Let c > be a sufficiently 
large constant, e A, Qi E M., qi > cqi+i, for i — 1, . . . , |.4|. Define lifting h : A ^ M^, where 
h{ai) = {wi, Qi), for i = 1, . . . , m, and h{ai) = (0, qi), for i = m + 1, . . . , \A\. Then, projecting the 
upper hull of A'^ to M^" yields the triangulation of A obtained by the 2-step construction. 

This is the implemented method; although different from the perturbation in the proof of 
Lem. [5] it is more efficient because of the reuse of triangulation Tq in Step 1 above. Moreover, our 
experiments show that it always validates the two conditions in Sect.|3] 

Fixing the dimension of the triangulation at compile time results in < 1% speedup. We 
also tested a filtered kernel with a similar time speedup. On the other hand, triangulation is 
expected to implement incremental high-dimensional regular triangulations wrt to a lifting, faster 
than the above method [TP] . Moreover, we use a modified version of triangulation in order to 
benefit from our hashing scheme. Therefore, all cells of the triangulated facets of 77 have the same 
normal vector and we use a structure (STL set) to maintain the set of unique normal vectors, 
thus computing only one regular triangulation per triangulated facet of U. 
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We perform an experimental analysis of our algorithm. We design experiments parameterized 
on: the total number of input points |^|, the dimension n of pointsets A^, and the dimension 
of projection m. First, we examine our algorithm on random inputs for implicitization and u- 
resultants, where m = n + 1, while varying We fix (5 G N and select random points on the 

(5-simplex to generate dense inputs, and points on the ((5/2)-cube to generate sparse inputs. For 
implicitization the projection coordinates correspond to point an = (0, . . . , 0) G For n — 2 the 
problem corresponds to implicitizing surfaces: when \A\ < 60, we compute the polytopes in < Isec 
(Fig. 4(a)). When computing the u-resultant polytope, the projection coordinates correspond to 

= {(1, . . . , 0), . . . , (0, . . . , 1)}. For n = 2, when \A\ < 500, we compute the polytopes in < Isec 
(Fig. [4^01. 

By using the hashing determinants scheme we gain a 18x speedup when n = 2,m = 3. For 
TO = 4 we gain a larger speedup; we computed in < 2min an instance where \A\ = 37 and would 
take > Ihr to compute otherwise. Thus, when the dimension and |^| becomes larger, this method 
allows our algorithm to compute instances of the problem that would be intractable otherwise, as 
shown for n = 3, TO = 4 (Fig. |4(b")] ). 

We confirm experimentally the output-sensitivity of our algorithm. First, our algorithm always 
computes vertices of U either to extend U or to legalize a facet. We experimentally show that 
our algorithm has a subexponential behaviour wrt to both input and output (Fig. 4(c)[ |4(d)[ ) and 
its output is subexponential wrt the input. 

We explore the limits of our implementation. By bounding runtime to < 2hr, we compute 
instances of 5-, 6-, 7-dimensional 77 with 35K, 23K, 500 vertices, resp. (Table [I]). 

We also compare with the implementation of [21], which is based on Gfain library. They 
develop two algorithms to compute projections of N(TZ). Assuming TZ defines a hypersurface, 
their methods computes a union of (possibly overlapping) cones, along with their multiplicities [211 
Thm.2.9]. From this intermediate result they construct the normal cones to the resultant vertices. 
We compare with the best timings of Gf an methods using the examples and timings of [5T]: our 
method is faster for m < 7, \A\ < 16 and competitive (up to 2 times slower) when to = 5, \A\ = 20 
or TO — \A\ = 15. We run extensive experiments for n = 3, considering implicitization, where 
TO = 4 and our method, with and without using hashing, is much faster than any of the two 
algorithms based on Gf an (Fig. |4(b")] ). However, for n = 4, to = 5 the beta version of Gf an we run 
was not stable and always crashes when |^| > 13. 

We analyze the computation of inner and outer approximations Q and ■ We test the 

variant of Sect.jsjby stopping it when ^^|^^^^^ > 0.9. In the experiments, the number of Q vertices 

is < 15% of the U vertices, thus the speedup is up to 25 times faster than the exact algorithm 

and ^l^^-* < 1.04 (Table pi. The bottleneck of this computation is the computation of vo\{Q^), 

which is in H-representation. We use Irs through polymake in every step to compute vo\{Q^) 
because we are lacking of an implementation that given a polytope P in H-representation, its 
volume and a halfspace H, computes the volume of the intersection of P and H. Note that we 
don't count this computation time in total time. 

Next, we study procedures that compute only the V-representation of Q. For this, we count how 

many random vectors uniformly distributed on the TO-sphere are needed to obtain ^p/^^j^j > 0.9. 
This procedure runs up to 10 times faster than the exact algorithm. Its drawback is that it does 

not provide guarantees for ^'^Y'^^ because we do not know U in advance. To overcome this, we 
^ ° vol(77) ' 

test a variant procedure which starts with a sequence of random vectors and produces vectors 
that are affine combinations of the existing ones. It stops when the produced vectors do not give 

a new vertex for Q. Even if it computes Q s.t. > 0.9 up to 10 times faster than the exact 

algorithm, on average it computes Q s.t. ^°|^^^^^ > 0.8 in similar runtimes (Table j2j). 
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6 Future work 



As shown, polymake's convex hull algorithm is competitive, thus one may use it for implementing 
our algorithm. On the other hand, triangulation is expected to include |10| fast enumeration 
of all regular triangulations for a given (non generic) lifting, in which case U may be extended by 
more than one (coplanar) vertices. 

The resultant edges are associated to flips on mixed cells [29] . We studied them algorithmically 
[l2] and may revisit them to generate all neighbors of a vertex of N{TZ). 

Any incremental convex hull algorithm has a worst-case super-polynomial total time complexity 
[5] in the number of input points and output facets. The question is whether there is a polyno- 
mial total time incremental algorithm for U. For this we study the structure of N{TZ) through 
experiments; it appears to be an exciting family of polytopes. 
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